The Greisen Equation Explained and Improved 
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Analytic description of the evolution of cosmic ray showers is dominated by the Greisen equation 
nearly five decades old. We present an alternative approach with several advantages. Among the 
new features are a prediction of the differential distribution, replacing Greisen's form which fails to 
be positive definite. Explicit comparison with Monte Carlo simulations shows excellent agreement 
after a few radiation lengths of development. We find a clear connection between Monte Carlo 
adjustment of Greisen's form and underlying physics, and present a concise derivation with all steps 
explicit. We also reconstruct the steps needed to reproduce Greisen's approximate formula, which 
appears not to have been published previously. 

PACS numbers: 96.50.sd 



I. INTRODUCTION 

The calculation of cosmic ray showers dates back to 
1937, when Carlson & Oppenheimer [l[ and Bhabha & 
Heitler |2| developed the theory of shower evolution. Un- 
til computers became powerful enough to run simula- 
tions, only analytic methods were available. The Greisen 
equation from 1956 i3J became the most popular solution. 
While state-of-the-art Monte Carlo simulations can yield 
very accurate results, analytic expressions are still used 
extensively. An analytic formula is invaluable for "back- 
of-the-envelope" calculations and order-of-magnitude es- 
timates. Analytic expressions can also be central to sim- 
ulations: by combining Monte Carlo calculations at high 
energies, where "new physics" happens, with analytic so- 
lutions for the lower energies, computing time can be re- 
duced enormously. 

Here we present the Greisen equation and show — as 
far as we know — its derivation for the first time in the 
literature. We then introduce an alternative treatment 
of the physics, the electro-photon model, which allows 
simple clear treatment without need for the ad-hoc re- 
placements used by Greisen. We show that the electro- 
photon (EP) equation agrees well with the Greisen equa- 
tion in its domain of applicability, and works much better 
where Greisen fails. Strictly speaking Greisen's formula 
applies to the integral distribution Tl(t, Eq/e c ), the num- 
ber of particles in the shower at depth t, given the energy 
of the incoming particle Eo and the critical energy e c . 
The corresponding differential distribution n(t, Eq, E) 
can also be obtained from an equation given by Greisen. 
Tr(t, Eq, E)dE is the number of particles with energies 
between E and E + dE. It turns out that for energies 
close to the energy of the incoming particle, the Greisen 
differential distribution is negative, as shown in detail 
below. Even though Greisen did not claim that his equa- 
tion is valid in the range where the number of particles 
is small, it is a disturbing fact that his formula has the 
feature of coming from an unphysical distribution. 

Yet Greisen's formula is simple and works amazingly 
well for many purposes. Modern Monte Carlo results are 



invariably presented via formulas marking up Greisen's 
form with effective parameters or variable substitutions. 
This suggests there is physics lurking in the empirical 
re-scalings. We find it by first deriving an improved for- 
mula, valid for both the differential and integral distri- 
butions of a shower, in which the electron- and photon- 
components are merged into one effective degree of free- 
dom. Naturally this opens up the opportunity to tune 
what is meant by an "electro-photon" with numerical 
work. Meanwhile the improved formula is as simple as 
Greisen's, while avoiding the unphysical elements and be- 
ing free of ad-hoc steps. For a cosmic ray with primary 
energy Eq, our practical formula is: 



U(t, E /e c ) = A EP e- tEP h (2\/^mT) -(2) 



Here /3o = \n(Eo/e c ), and Iq, I\ are modified Bessel 
functions. The single-species effective depth tEP is given 
by 



t EP = Xt + a EP + [3 Q (l - A) - 1/2 



(3) 



where A = 2/V3 « 1.155, parameters a E p and the 
normalization Aep are given in Table H These Monte 
Carlo-based re-scalings are perfectly consistent with the 
model's formulation as due to merging two species of par- 
ticles into one. 

Readers concerned mainly with practical results can 
compute with Eqs. [TJ [21 which produce rather fine agree- 
ment with Monte Carlo simulations. In the rest of the 
paper we show how the results were obtained. Section 

II recounts the reconstruction of Greisen's steps. Section 

III shows how the well-behaved analytic formula is ob- 
tained. Section IV compares the electro-photon model 
to Monte Carlo results, explains the need for re-scalings, 
and contains Table |U 
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II. THE GREISEN EQUATION 

Greisen's approximation for the number of particles 
IT (t, Eq/ec) in an electromagnetic shower "in the region 
where the number of particles is large" [|J , long referred 
to as the "Greisen equation" , is: 



U(t,E /e c )^^e< 1 -^). 



(4) 



Here t is the distance the shower has developed in units 
of the radiation length, /3q is given by 



00 = In 



E 



(5) 



where Eq is the primary energy, e c is the critical energy. 
The variable s came to be called the shower age, given 
by 



3/ 



t + 2(3 



(6) 



A. Derivation of the Greisen Equation 

Unfortunately, Greisen states this equation without 
giving a derivation. Most of the papers and textbooks 
that use the Greisen equation cite Greisen's 1956 paper, 
without giving the full derivation either (see, e.g. [4j, |5[). 
We have reverse-engineered the steps that Greisen took 
to find his equation and present them in this section. 



1. Rossi and Greisen, 1941 

A good starting point for the Greisen equation turns 
out to be the article by Bruno Rossi and Kenneth 
Greisen Q (RG). It summarizes the state of cosmic 
ray shower physics at that time. The important part 
for the Greisen equation is the "Approximation B" for 
shower evolution. Approximation B includes pair pro- 
duction, bremsstrahlung and collision energy losses of 
the electrons, but neglects the Compton effect. RG start 
with the asymptotic formulae for Bremsstrahlung and 
pair production, 



dcr 1 

dk ^ k 



dcr 
dE 



k 
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l 



x(l 



(7a) 
(7b) 



where k is the photon energy, E the electron energy, 
y = k/E and x — E/k. Using Eqs. [7] and an expression 
for collision energy losses, RG set up the shower evolu- 
tion equations. They can be solved in a formal sense an- 
alytically, using Mellin integrals. We review this shortly: 
After several pages of calculations, a certain saddle point 



approximation and dropping a negligible term, the result 
is 



n(t) = 



1 H^K^-s) fE 



2tts y/\'{{s)t+ l/s 2 



±2 p \As)t 



where 



t = - 



\[(s) 



In 



Er 



(8) 



(9) 



-ffi(s), Ki(s, — s) are rather complicated functions that 
are tabulated in the article and Ai(s) is an analytic func- 
tion, also given in the article. 



2. Greisen's approximations 

Several approximations must be done to get from RG 
to the Greisen equation. The primary idea is to focus on 
the behavior of the shower around the shower maximum, 
that is, for shower age parameter s = 1. 

First, the function Ai(s) is approximated. A Taylor 
expansion of sA^(s) around s — 1 yields: 

sAi(s) « -0.938 + 0.460(s - 1) - 0.037(s-l) 2 + 

+0((s-lf). (10) 

By reasonable guesswork we are able to deduce that 
Greisen replaced coefficients in the following way: 

0.938 -► 1 

0.460 -► 1/2 

0.037 -► 0. 

This yields: 



X 1 (s) 



(s- 1 -31ns). 



(11) 



In Eq. ©, we can now drop the l/s term (since 0q 3> 
l/s). Using Eq. (fTl"j) we then get 



3f 



(12) 



which is just Eq. ([6]) used in the Greisen Equation. 

With the help of both Eq. (JTTJ) and Eq. O and a 
little bit of algebra we can rewrite a part of the Rossi 
and Greisen equation, Eq. ©: 



Eq 

Sr. 



,Ai(a)i _ i(l-flna) 



(13) 



In the remaining part of Eq. ([5]), we drop the l/s 2 term 
since it is small compared to X'{(s)t, and evaluate the rest 
at s = 1: 



1 H^Kfa-a) 



2ns 



0.3162 0.31 



(14) 
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Putting eqs. (| 13|) and fl4]) together, we get the Greisen 
equation. 

We note that in some places s is evaluated at s = 
1 whereas at other places it is left as a variable. The 
replacements above are not controlled approximations: 
the Mellin variable s is supposed to be integrated over 
a certain contour, and there is no way to predict the 
effects of the substitutions made. The particular steps 
above are definitely not designed to work in regions where 
the distribution deviates from its maximum. Notice also 
that the method of Mellin moments applies to both the 
differential distribution and the integral particle number, 
but the steps used are not equally good for both. These 
facts begin to explain the problem mentioned earlier that 
the differential distribution corresponding to Greisen's 
formula is unreliable. 



n 








FIG. 1: Differential distribution as obtained from the Greisen 
equation before, at and after shower maximum (arbitrary 
units) ; energy of the primary is 1 x 10 9 . 



B. Differential Distribution 

The Greisen equation does not provide a differential 
distribution, but in Greisen's article [3J], we find the 
following: "One may compute the approximate num- 
ber of particles having an energy exceeding E, provided 
Wq 3> E Eq (Eq E 3> e c in our convention) 
from Eq. 1 (the Greisen equation) by replacing the coef- 
ficient 0.31 with 0.135 and substituting for f3 the quan- 
tity (3 = ln(W Q /E) (f3 = ln(E /E) here): 7 It might seem 
that these substitutions come from nowhere and are com- 
pletely arbitrary, but it can be shown that they follow 
from Rossi and Greisen, 1941 [(|, by just the same steps 
as the Greisen equation. The new equation reads: 



n (t, E' > E) 



0.135 



t(i- 



'). 



(15) 



This allows us to calculate the differential distribution 
7r(i, E), that is the number of particles with energy E at 
shower depth t: 



n(t, E) 



dll(t, E' > E) 
dE ' 



(16) 



It turns out Greisen's differential distribution is negative 
for energies close to the energy of the incoming parti- 
cle, as shown in Fig. [1] Needless to say, this behavior of 
7r(t, E) is unphysical. Even though the equation is not 
supposed to be used at energies close to the energy of 
the incoming particle, this unphysical behavior is a dis- 
turbing fact of the replacements and substitutions made. 



C. Greisen and Monte Carlo 

Comparisons of the Greisen Equation with Monte 
Carlo calculations are very common. In order to make the 
Greisen equation agree with the simulations, Fenyves et 
al. Q introduced two parameters: a(E) shifts the Greisen 



equation along the radiation length axis, and A{E) re- 
scales the Greisen equation by a factor. Then 



n(t) 



Q.ilA{E) 



,ti(l-f lnsi) 



(17) 



where t x = t + a{E) and si = 3<i/(ti + 2/3 ). 

It is impossible to summarize all comparisons of the 
Greisen equation with the Monte Carlo calculations here. 
One of the more insightful results was obtained by Sci- 
ascio, Piazzoli and Iacovacci {SPI) SPI fit a(E) 
and A(E) for different threshold energies for both the 
electron component and photon component of electro- 
magnetic showers. SPPs paper then finds that the data 
points from the Monte Carlo lie exactly on the lines from 
the (adjusted) Greisen equation. Given the number of 
unjustified steps and substitutions this is indeed a puz- 
zle. 



III. THE ELECTRO-PHOTON APPROACH 

In the high energy limit all particles act as if mass- 
less and the distinction of electron-positron pairs versus 
photons ceases to be physically meaningful. This sug- 
gests we should drop the distinction in the mathematics. 
The basic assumption for the electro-photon approach is 
to consider only one species of effective particles, which 
for discussion we designate the massless electro-photons. 
They will replace the electrons, positrons and photons in 
the regular shower models. However the limiting cross 
sections of pair production and bremsstrahlung arc not 
quite equal. With one species replacing two, what is 
meant by the "radiation length" must be adjusted to a 
suitable value. We gain the freedom to fit the scaling 
parameter for the effective radiation length after solving 
the model. 
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FIG. 2: Energy conservation in the SP-model (the radius 
represents the energy of the particles): an incoming EP of 
energy E (A) loses energy by radiating a "bremsstrahlung" 
EP of energy E 1 (B); A' should now have energy E — E' (solid 
line), but in the _EP-model we end up with 1 — E / ' E electro- 
photons of energy E (dashed line); obviously, this conserves 
energy. 



A. Setting up the equations 

There is little freedom in the cross-section for EPs due 
to dimensional analysis of the high energy limit. Quan- 
tum electrodynamics is scale free, so that: 



da 1 
dk^k' 



(18) 



This is also an approximation to the Bremsstrahlung 
cross-section, Eq. ([7a]) . with 



;V + V 



1. 



(19) 



Using this cross-section we can write down an evolution 
equation for the EPs, 



d7r(i, E) 
~dt 



f°° 1 
J dE w(t, E) — 



(20) 



where t is the shower depth in radiation lengths. 

This, of course, describes only the gains and not the 
losses and violates energy conservation. Subtract 7r(i, E) 
on the right hand side to get 



dir(t, E) 
dt 



,oo 1 

J e dE'7r(t, E')--n(t,E), 



(21) 



which incorporates energy conservation. The facts of ex- 
act energy conservation will be more easily seen once we 
have the Mellin transform of this equation. 

The loss part of Eq. [2D can be pictured as shown in 
Fig. [2j assume one incoming electro-photon of energy E. 
It radiates an EP of energy E' . Then the original EP 
should have only energy E — E' left. However, according 
to the above equation we end up with 1 — E'/E electro- 
photons with energy E. Obviously, this conserves energy. 
And since in most cases E'/E <C 1 (due to the 1/E 
behavior of the cross-section), this approximation is quite 
reasonable. 



B. Solution of the electro-photon model 



To solve the evolution equation, Eq. (|2Tj) . we take its 
Mellin transform (see Appendix [A"| . The only non-trivial 
part is the integral on the right hand side: 



rdE N r dE ,]_ n{E/) 

Jo E J E E K ' 

poo pE' 

= / dE'n(E') / dEE 
Jo Jo 



iN-1 



N - 1 



w(N) 



(22) 



where 7r(A^) is the Mellin transform of k{E). Therefore 
the evolution equation is: 



d7r(t, N) 



1 



- 1 7r(t, N). (23) 



dt \N-1 
The N — 2 mode gives the total energy, and we have 
dn{t, 2) ( 1 



df 



2- 1 



1 7r(t, 2) - 0. (24) 



Therefore energy is conserved, just as expected. 

The evolution equation can be easily solved by an ex- 
ponential 



7r(*, N) = e(^T" 1 ) t 7r(0, N) 



(25) 



where 7r(0, N) describes the incoming particle. We as- 
sume one incoming electro-photon with energy Eq: 



and 



n(t = 0,E) = 8{E - E a ) 



7r(0, N) = J™ ^§E N S(E Eq) = E^-\ 



(26) 



(27) 



Therefore, the Mellin transform of the EP distribution 
is 



7r(t, N) = ei^-^E^- 1 . 



(28) 



Now we just have to transform this back to 7r(i, E). 
Taking the inverse Mellin transform yields: 



1 



a-Hoo 



nit, E) = / dNE- N K(t, N) 

2ki 



1 

2kI 



a — ioo 



dJV J B- iV e(^- 1 )*^- 1 



a — i oo 



2niE 



dNe^^ N -^ ln ^. (29) 



For \n(Eo/E) > 0, we can close the contour in the left 
half-plane and get 

n{ t 7 E ) = ^5Hes ( e *£i+("-i)to^ . (30) 
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FIG. 3: Differential distribution for the electro-photon model 
before, at and after shower maximum (arbitrary units); en- 
ergy of the primary electro-photon is 1 x 10 9 . 



FIG. 4: Dependence of differential distribution on shower 
depth in the electro-photon model for certain fractions of the 
energy of the incoming particle. 



The residue must be evaluated at the essential singularity 
N = 1. It is the residue of an exponential of the form 
exp(a/a; + fix) which is treated in Appendix IbI Using 
the results, we get 



7t(t, E) 




(31) 



where I\ is a modified Bessel function. 

For ln(E /E) = 0, that is E a = E, we cannot close 
the contour in the left half-plane 9]. Instead, we go back 
to the evolution equation, Eq. 1211 and using the initial 
conditions, Eq. [20 we see immediately that 



ir(t, E = E ) = e - t S{E - E ). 
Together, we get 



(32) 



7r(t, E) 



iln§ 



for E < E 
for E = E Q , 



In 

^ (33) 

The exponentially decaying 5-function for E = Eq can 
be understood as the probability for the original particle 
not to interact down to a certain radiation length. How- 
ever integrating the ^-function amounts to less than one 
particle, which compared to the large number of particles 
in the shower is negligible. Therefore we will drop this 
term from now on. 

In Fig. [3] we show that the differential distribution is 
everywhere non-negative as expected. Fig. U shows the 
i-evolution for the EP distribution function for certain 
fractions of the energy of the incoming particle. 

To get the total number of electro-photons in the 
shower, we integrate the differential distribution with re- 
spect to E. In principle a calculation includes collision 
losses: particles lose the critical energy e c per radiation 
length due to collisions. Therefore, particles with en- 
ergies less than e c will drop out of the shower rapidly. 
Hence, we will only count particles with energies greater 



/ \ 





Eo 


= 10 19 eV 




Eo 


= 10 18 -' eV 




Eo 


= 10 18 eV 




FIG. 5: Dependence of total particles in the shower on shower 
depth in the electro-photon model for E = 10 19 eV, 10 18 5 eV 
and 10 18 eV; the critical energy is e c = 81 MeV. 



than e c : 

n(t) 



&En{E,t) = e-'/o f 2Wtln^ j . (34) 



Here Iq is a modified Bessel function and e c is the critical 
energy. A plot of H(t) for three values of the energy of 
the incoming EP is shown in Fig. [5J 

There remains the task of assigning an effective depth 
t — * t ep and making use of the numerical results to cor- 
rectly scale our equations. 



IV. ELECTRO-PHOTON MODEL AND MONTE 
CARLO SIMULATIONS 

We now compare the results from the electro-photon 
model to Monte Carlo simulations. To do so, three pa- 
rameters will be fit. One is an overall scaling factor and 
the second a shift in i-direction. These two steps have 
been performed previously to make the Greisen equation 
fit the Monte Carlo simulations, so they need no further 
justification. 
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Earlier we mentioned that the effective cross sections 
of EPs must be a compromise between pair-production 
and bremsstrahlung processes. We therefore complete 
the EP model by adjusting a single effective radiation 
length, present from the start, as a free parameter. 

We want to make use of the existing fit parameters 
that match the Greisen equation with Monte Carlo re- 
sults. To do so, we will find the parameters that connect 
the electro-photon model with the Greisen equation, us- 
ing analytic methods. Since for most applications the 
energy of the primary will be much greater than the crit- 
ical energy, the argument of the Bessel functions will be 
large and we can use the asymptotic formulae 



h(x) 




(35) 



First, we require that the maximum occurs at the same 
shower depth. For the Greisen equation, the shower 
maximum occurs for i max = /3o- Using Eqs. [35l we get 
tmax — Pa — 1/2 for the EP-model. Second, we require 
that the number of particles be fixed at shower maxi- 
mum. For this it is sufficient to keep only the leading 
term in Eqs. 1351 This second requirement leads to the 
scaling factor of A = 0.31 x y/An sa 1.099. Third, we 
adjust the effective radiation length. This is done by 
requiring that the second derivatives of the number of 
particles with respect to shower depth agree at shower 
maximum. This yields a factor for the effective radiation 
length of A = 2/^3 » 1.155. 

Incorporating these results, we get the parametrized 
EP equations: 



7T(i, E) = 



A E pe~ 



tEP 



't EP \n 



Eq 



(36) 



E V In % \ V E 
II(i, E /s c ) = A EP e- tEP I (2VW0) • (37) 
Here the adjusted shower depth is given by 

t E p = \t + a E p + - X) - 1/2. (38) 

A E p and a E p can be obtained from the parameters A(E) 
and a(E) in the Fenyves et al. Q parametrization of the 
Greisen equation in the following way: 



A EP = AA(E) 
a E p = Xa(E) 



(39) 

where A = 0.31 x y/I^ w 1.099 and A = 2/V3 « 1.155 as 
above. 

With the numerical data obtained by Sciascio, Piazzoli 
and Iacovacci @, we calculate the parameters A E p and 
a E p for the electro-photon model that agree with these 
simulations. The parameters are summarized in Table [U 

Obviously, when fitting new Monte Carlo simulations, 
it is not necessary to take the detour via the Greisen 



1 3 TcV 




FIG. 6: Dependence of number of particles in the shower 
on shower depth for E = 10 3 TeV and 10 2 TeV for electro- 
photon model (lines) and Monte Carlo [H (squares and tri- 
angles). Critical energy e c — 81 MeV; threshold energy 
E th = lMeV. 
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FIG. 7: Dependence of number of particles in the shower 
on shower depth for E = 10 3 TeV, 10 2 TeV, 10 TeV, 1 TeV 
and 0.1 TeV for electro-photon model (lines) and Monte Carlo 
[i| (crosses). Critical energy e c = 81 MeV; threshold energy 
E th = lMeV. 



equation and A E p and a E p can be directly fit from the 
data. 

Incorporating the parameters into the electro-photon 
model, we get the results cited in the Introduction, and 
shown in Figures [5] and \7\ Looking at the plot with the 
linear scale, Fig. [51 there is beautiful agreement between 
the electro-photon model and the Monte Carlo simula- 
tions. The plot with the logarithmic scale, Fig. [3 shows 
that development over a few radiation lengths is needed 
to converge to the simulations, just as expected. The 
agreement is very good well before the shower maximum 
and right across the tail of the shower. 



V. CONCLUSIONS 

Greisen's 1956 formula cannot really be derived from 
coherent mathematical approximations, and it must have 
been motivated from early numerical work in a time pre- 
dating high-speed computers. The underlying physics is 
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TABLE I: Aep and asp for the electron component and pho- 
ton component of showers, depending on the threshold energy 
E t h- Values are based on data obtained by Sciascio, Piazzoli 
and Iacovacci 



E th (UeV) 


electron component 

Aep Qjep 


photon component 
Aep clef 


1 


1.01 


0.00 


5.27 


-1.02 


5 


0.82 


0.22 


3.27 


-0.80 


10 


0.69 


0.40 


2.34 


-0.66 


15 


0.59 


0.52 


1.88 


-0.52 


20 


0.55 


0.61 


1.59 


-0.42 


50 


0.35 


0.96 


0.81 


0.14 


100 


0.23 


1.39 


0.45 


0.73 



extremely simple, and dominated by the 1/k dependence 
of the cross sections, which is a consequence of the scale 
invariance of quantum electrodynamics. A treatment 
retaining nothing but the leading 1/k dependence and 
combining electrons and photons into one common entity 
leads to simple analytic formulas with many advantages. 
First, the derivation is concise and explicit, and there 
is no need to replace expressions by proxies. Second, a 
usable differential distribution is obtained, which is nec- 
essary for many purposes. Finally the results broadly 
agree with Monte Carlo simulations provided that an ef- 
fective radiation length is employed, as consistent with 
the model. Corrections to this starting point (none are 
indicated numerically) would probably involve the com- 
plicated functions of ratios of two scale parameters seen 
in the exact solutions and Approximation B. 
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APPENDIX A: MELLIN TRANSFORMS 

1. Definition 

The Mellin transform f(N) of a function f(E) is de- 
fined by 



/(AO 



(Al) 



2. Inverse Mellin Transform 

The inverse of the Mellin transform is given by 

-i />a+ioo 

f(E) = — / dNE- N f(N) 

= -L f dNE- N f(N). 
2m J c 



(A2) 



The contour of integration C is parallel to the imaginary 
axis. 
Proof: 



dNE~ N f(N) 
- [ dNE- N r^L E > N f(E') 



1 

2 

- SiffW^™ 



Wc know that 



— / dNe Nh ^ E '^ E ^ 

27Ti Ja-ioo 

= ^-e aX < E '' E h ( +0 ° dxe lxl ^ E ''^ 

27T« J-oo 



(A4) 



and with this, we have: 



= f(E). (A5) 
APPENDIX B: RESIDUE OF exp(a/x + f3x) 



To calculate the residue of 

ei +f>x 

we write the exponential as a series: 



(Bl) 



'-+/3x _ 



2 = 

i=0 j=0 XJ/ 



(B2) 



For the residue we only need the values of j such that 
-3 + j) = -1, i-e- j = {i + l)/2. So let i = 2k + 1 
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which yields j = k + 1 . Then where 1\ is a modified Bessel function. 
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